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ABSTRACT 

In this paper we propose a new method for reconstructing the surface density of matter in 
flat disks of spiral galaxies. The surface density is expressed through observational rotation 
velocity curves of visible matter in the disks of spiral galaxies. The new method is not based 
on quadrature of special functions. The found solution is used for processing and analysis 
of observational data from several spiral galaxies. The new method can be used to more 
accurately estimate the amount of dark matter in spiral galaxies. 

I. INTRODUCTION 

Let us consider a model of a spiral galaxy with all matter distributed in a thin disk (thickness 
disk is neglected). Accounting for the components of a spherical halo of the galaxy see in section IVTl 

Let us consider the axisymmetric model in cylindrical coordinates {p,z,(p). 

In the disk we assume a stationary rotation of dust-like matter, i.e. equality (modulo) of the 
gravitational and centrifugal forces: 

dp p 

Here V{p) is the rotation velocity curve of the matter and ^{p,z) is the Newtonian gravitational 
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potential of the dust matter, which may be expressed as: 

oo 7r 

$(p,z) = -G [ RdR [ dip ^ ''^^^ = (2) 

J J \ +E? - 2pRcos in + 







where cr{p) is the distribution of the surface density of matter in the disk. 

We pose the problem: Express the distribution of the matter surface density, cr{p), in the disk 
from the rotation velocity curve of the matter V{p)- 

The potential ^{p,z) is a function of cr{p). Gauss's theorem for this potential leads to the 
expression: 

Hp,z),,\^^^ = 27rGa{p) (3) 

Here G is the gravitational constant. 

Thus, we can obtain an integral equation expressing the velocity V by the integral of the surface 
density a. This equation can be inverted (see to express the surface density a through the 

rotation velocity curve in the disc: 

oo oo 

a{p) = J dkJo{kp)a{k), a{k) = ^J dR Ji{kR)V^{R) (4) 



Here Joix) and Ji(x) are Bessel functions. 

Using this method to determine the surface density cr{p) has significant shortcomings: 

1. To find the desired function cr{p) it is necessary to integrate a fourfold integral, whereas the 
Bessel functions themselves are each expressed in terms of integrals. 

2. In the integrand of Q are oscillating and slowly decaying (as cosx/^/x) Bessel function, 
which also makes it difficult to integrate. 

3. The integrals has infinite upper limits and therefore they requires a large numerical resources 
(change of variables leads to new difficulties due to the oscillating integrands). 

All of these deficiencies substantially affect the quality of the results (obtained by numerical 
methods) . 



In papers 



a more convenient method based on integrals of elliptic functions. It was shown 



that the main errors arise due to ignorance of the distribution rotation velocity curve V{p) at 



arge 



values of p. The same problem (stability of the numerical calculation of a) was studied in 

In this paper we propose a new method for solving the problem of determining the surface 
density cr{p) from the disc velocity curve V{p). Our method does not require knowledge of the 
numerical values of special functions (Bessel functions and/or elliptic functions). 
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II. THE NEW METHOD 

Using the relation Ji(x) = —Jo{x),x for the Bessel functions and applying integration by parts 
to the integral for (^{k) we rewrite expression 

oo oo 

2TrGa{p) = J dk J dR Jo{kp)Jo{kR) [V^{R)],n (5) 


Here we took into account that 



oo 



/I 1 If 

dR -Jo{kR),R V\R) = - -Jo{kR) V\R) + ^ dR ,h{kR) [V^{R)U 

° 

Also we took into account that at x = the function Jq{x) is limited and V'^{x) vanishes (no 
velocity on the axis). At x — )• oo the opposite is the case, i.e. the function ^■^(x) is limited and the 



function Jo{x) decays as cos(x — 7r/4)y^2/(7rx). 

We denote the gradient of the square of the velocity as f{R): 

f{R)^[V\R)U (6) 
Substituting the integral definition of Bessel functions of zero order 

Jo(x) = iy"e^^''=°^"da (7) 



into eq. (0), gives: 

7r 7r oo oo 

2TxGa{p) = ^ j da j dp j dk j dR e^kpcosa^ikRcos/s j^j^-^ 

-oo 

It was here taken into account that Jo{x) is an even dunction, so that the integration over k in 
([5]) can be made by changing the limits of integration to go from — oo to oo and multiplying the 
result with a factor ^ . 

Let us now consider the integral representation of the Dirac delta-function: 

oo 

6{x) = — [ e'^'^dk 
2vr J 

— oo 

and substitute it in (l8|): 

7r 7r oo 

2TTGa{p) = ^ j da j d(3 j dR 6{pcos a + Rcos p) f{R) (9) 
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Here we take into account the integration rules of the complex delta-function: 

\9{R),R I 

where Rq are the simple roots of the equation g{Ro) = 0. Also we should have that Rq > 0, which 
means that the integration region with respect to a and /3 reduces to the two sectors: 

f a G [0; 7T/2] ^ (ae [tt/I; vr] ^^^^ 
■ \ /3 e [7r/2; tt] ' \ /? G [0; 7r/2] 

By renaming the variables a and /3 in these sectors^ we may reduce the integral ([9]) to the form: 

2.Gaip) = Uda7dpl^, Ro = p'-^ (11) 
TT J J smp smp 


This expression, (jlip . is the desired result which is without the deficiencies mentioned in section HI 
Note that in expression (jlip we used the so-called holographic principle, i.e. each point of the 
density distribution cr{p) contains the information about all points of distribution f{R). This is 
due to integration over the angles a and /3. Integration over the angle a provides a scanning range 
[0, p], and integration over the angle /3 provides a scanning range [p, oo). In addition, sin/3 in the 
denominator of the integrand provides increased weight of the integrand f{R) at large R. The 
converse is also true; each point of the distribution f{R) contains information about all points of 
the density distribution a{p) - see section Hill 

This interesting feature is characteristic for the mass distribution in a flat disk^. 
Using the formula V'^{Ro),p = V'^{Ro),Rq ■Ro,p = f{Ro)cosa/ sin/3 we can rewrite eq. (fTTj) in 
the form: 

n/2 7T/2 ^ p 

27rGa{p) = , U{p) = -[ da [ dp , u{p) = [ 27rGa{p) dp. (12) 

dp J J cos a J 



Using this expression it is convenient to obtain an expression for the total mass of matter inside 

the disk radius p by integrating eq. (jl2p by parts: 

p p 
GM{p) = j 2TrGpa{p) dp = U{p) ■ p- j U{p) dp (13) 



^ In the first case, we redefine a = a, /3 = 7r/2 + /3; in tiie second case: a = tx ~ a, P — 7r/2 — /3. 

^ In tiie case of spfierical symmetry, tiie model will no longer possess such a property; each point R in the distribution 

of velocity will carry information only about the total mass inside a sphere with a radius R, and value of the mass 

density will be associated only with the pressure of matter (not with velocity). 



5 




III. CHECKING THE METHOD 



For basic testing of our new method, we use the fohowing approach; We integrate the expres- 
sion ([T]) with respect to p: 

H,{p)^j^-^dp (14) 



On the other hand, this same expression can be obtained from ([5]): 

oo 7r/2 

H2{p) = C/(oo) - - fdR- 27TGa{R) == , 7 = (^/2 (15) 

J J A/(ii-/9)2 + 4i?psm^7 

Or as: 

00 7r/2 

H2(p) = - I dR- 27:Ga(R) [ dj 1 1 ^ 

vr y ^ ' ' \ ^{R-pf + ARp^m' 

The integrals psp and p6p do not have non-integrable singularities and can therefore easily be 
found numerically. 




6 



Expressions for the functions Hi{p) and i?2(p) ought to be identical, as they must be equal to 
the expression ^{p, z = — $(0, z = 0). Thus, we may judge the practical precision of our new 
method by evaluating and comparing the expressions for Hi{p) and H2{p) (using our new method 
to calculate cr{R)). 

To evaluate and compare these expressions, and thus to put our new method to use, we have 
written a simple numerical code, the details of which are described in Appendix A. In Fig. [1] is 
shown an example of such a comparison, where we have used the sample velocity function : 



l + (10x)3 



as the basis for data. In the figure this velocity function and cr{x) are calculated by using our new 
method. Also, of course, is plotted the two functions Hi{x) and H2{x) that we wish to compare. 

As can be seen from the figure (and as has been confirmed by detailed analysis of the figure) , the 
two functions Hi{x) and H2{x) agree very closely. This is also the case for other velocity functions 
we have tested. Thus we conclude that our new method is both valid and accurate, even when 
doing the integrations numerically. 



IV. SPECIAL CASES AND ASYMPTOTIC DISTRIBUTION OF THE DENSITY OF 

MATTER 



In some special cases, expression (llip can be reduced to an ordinary quadrature and we may 
find the asymptotic behaviour of cr{p). 

1. Consider the simplest case for the distribution function f{R) and the rate V{R): 



fiR) 



V{R) 



VATT? for R e [0; Ri] 
V/o • ^1 for i? G (i?i; oo) 



(18) 



/o for i? G [0; 
for i? G (fli; oo) 

The integral over f3 in expression (jlip can be easily evaluated and we may obtain the quadrature: 

7r/2 

da, ai = < 



2TTGa{p) = ^ / In 
vr 



ai 



1 + a/1 - {pcosa/RiY 



or in another form: 



2-KGa{p) 



2-Ri/o 



pcosa/Ri 



In (1 + Vl-x2)/a 



for p < Ri 
arccos(i?i//9) for p > i?i 



(19) 



T^P J y/1 - {xRl/pf 



dx , xi 



p/Ri for p < Ri 
1 for p > Ri 



(20) 



From expression (fT9]) we find that the asymptotic behaiour at p/Ri — )■ goes as cr(p) oc ln(i?i/p). 
From expression (pO|) we find that the asymptotic behaviour at p/i?i — >• oo goes as a{p) oc (Ri/p). 
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It is noted that the asymptotic behaviour of has a weak singularity at the origin. This 
is connected with the non-rigid rotation of matter around the central axis (for the case of rigid 
rotation, we should have V{R) oc at i? — 0) - see eq. (fTSll . 

2. We now consider another simple case, that of bounded rigid-body rotation: 



fiR) 



V{R) 



fo ■ R/Ri for R G [0; Ri] 
for i? G (i?i; oo) 

In this case the analogue of the integral (jl9p has the form: 

tt/2 

27rGa{p) = / \/l — (pcos a/ RiY da , a2 = 

71" J 



^■R for Re [0; Ri] 



0.2 



^ for i? G (i?i; oo) 



for p < Ri 
arccos(i?i/p) for p > Ri 



(21) 



(22) 



This integral is expressed through the elliptic integral of the 2nd kind, here the asymptotic behaviour 
at p/Ri — )• is nonsingular and we have that cr{p) — )• const. At infinity, the asymptotic behaviour 
remains the same: That is, at p/Ri — t- oo we have cr{p) oc {Ri/p). 



V. OBSERVATIONAL DATA 



By using the method described above we have restored the matter surface density profiles for the 
discs of 4 spiral galaxies basing on their observed rotation curves. The four galaxies are : a) NGC 
2841, b) NGC 7217, c) NGC 7331, and d) NGC 5533 (the general characteristics of the galaxies 
are listed in Table [J and the restored surface density distributions are presented in Fig. [2|). 

We use for computing the two limiting cases (see Fig. [2|) for the behavior of the curve V{R) at 
large R > Rmax (where the velocity distribution is unknown): 

1. V = V{Rmax) = const 

2. V = V{Rmax) • \/R max 

/R 

The second case corresponds to a Keplerian tail in the velocity distribution of matter, quantities 
corresponding to this case, marked on Fig. [2] by a prime. 

The galaxies under consideration are all giant spiral galaxies of the Sab-Sb type and all of 
them are classified as unbarred This is important because the axisymmetry of the galaxies signifies 
regular circular rotation of the gas, providing a possibility to measure their rotation curves from a 
single spectral observation per galaxy if the spectrograph long slit is aligned with the disk isophote 
major axis. About accounting for the components of the spherical halo see section IVll The rotation 
curves have been obtained by measuring the line-of-sight velocities of the gas. Because the gas is a 
cold dynamical subsystem and it is confined to the thin disks, its rotation traces the gravitational 



FIG. 2: Observed and calculated data for the galaxies: 
(a)7VGC2841, {b)NGC7217, (c) NGC7S31 and (d) NGC5533. 

The horizontal axis is the distance in units of a; (x = R/ Rmax)- On the vertical axis plots: 

A) - V{x)/V{l); 

B) - M{x)/Mmax; C) - (T{x)/amax - for V {R > Rmax) = V (Rmax) = COnst; 
B') - M{x)/Mmax] C) - a{x)/a„,ax - for V{R > R^ax) = V {R^ax) ^/ Rmax I R\ 

D) — vertical line marking the end of the stellar disk to the 25-th blue isophotes. 
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(a) NGC2841: ij^a^i = 63 240Pc, 
A) V(l) ^ 294fcm/sec, 
B) M™„, = M(l) 1.22 ■ lO'^Me^ 
C) cr„„, = a(0) f« 3 432M0/Pc2, 
C") cr(0) ^ 3 406Mq/Pc2. 
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(c) NGC7331: R^a^ = 37 230Pc, 
A) V(l) ^ 238fcm/sec, 
B) M™, = A/(l) 4.93 ■ 1O"M0, 
C) o^a^ = a(0) f« 1883M0 /Pc^ 
C") cr(0) ^ 1854M0/Pc2. 




(b) NGC7217: R^a^ = 9 600Pc, 
A) V(l) 305fcm/sec, 
B) M„,, = M(l) ^ 1.94 • 1O"M0, 
C) o-„„, = o-(O) ^ 6 22OM0/Pc^ 
C") o-(O) 6O5OM0/Pc2. 
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(d) NGC5533: P™a^ = 77 875.2Pc, 

A) V(l) 22&km/8ec, 
B) M^a^ = A/(l) ^ 9.34 • 1O"M0, 
C) a^a^ = o-(O) ^ 2 877M0/Pc2, 
C") (j(0) 2 867M0/Pc2. 
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TABLE I: Global parameters of the galaxies 



NGC 


2841 


5533 


7217 


7331 


Type (NEDi) 


SA(r)b SA(rs)ab (R)SA(r)ab 


SA(s)b 


i?25, arcsec (LEDA^) 


244 


93 


117 


314 


i?25, kpc 


16.6 


19.3 


9.3 


23 


Distance, Mpc (NED) 


14 


43 


16.5 


15 


Mb (LEDA) 


-20.84 


-21.6 


-20.48 


-21.56 


Kotsini, km • (LEDA, HI) 


295 


206 


143 


237 


Mhi, (1O9M0)3 


2.8 


30 


0.7 


8.2 



^NASA/IPAC Extragalactic Database 
^Lyon-Meudon Extragalactic Database 

3from m for NGC 2841 and NGC 7331 and from M for NGC 5533 and NGC 7217 



potential distribution very well. For the inner parts of the rotation curves we have used spectral 
observations of the ionized (warm) gas which have a good spatial resolution (typically about 1 
arcec that corresponds to a few dozen of parsecs for 3 of the 4 galaxies). For the outer parts, we 
used radio observations of the neutral hydrogen at 21 cm which has a spatial resolution that is 
worse by an order of magnitude compared to the optical. However, the neutral hydrogen disks 
are usually more extended than the stellar ones which allows us to trace the rotation curves (and 
thus the gravitational potential distribution associated with them) much farther from the galactic 
centres than where the stars are seen. The literature sources which we have used to construct the 
combined rotation curves for our 4 galaxies are listed in Table HH 

It is interesting to compare the restored density distributions of the gravitating matter with the 
distributions of the stellar component in the galaxies considered; in all of them the stars are the 
dominating baryonic content. For this purpose we have used surface photometry data from the 
literature for our galaxies. NGC 7331 and NGC 2841 have extended radial profiles of the surface 
brightness in the NIR (2 mkm) photometric band K published in the survey by [lo| . the other 
two galaxies have been studied by us earlier in the ii-band (NGC 5533, see [111 ]) and in the I- 
band (NGC 7217, see jl2]). To derive a density distribution from the brightness distribution or to 
calculate the brightness distribution corresponding to some density distribution, we need to know 
the mass-to-light ratio in the specified photometric band appropriate for the stellar population of 
the particular age and metallicity. At the same time, the latter characteristics are only known 
for the disks of a few nearby spiral galaxies. For NGC 7217 and NGC 5533 we have studied the 
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TABLE II: References for the rotation curves and photometric profiles 



NGC 


References 


NGC 


References 


2841 


rot(Ha+[NII]): [16] 


7217 


rot(Ha+[NII]): [17] 




rot (HI): [18,i9j 




rot (HI): [9] 




density (HI): f8] 




density (HI): [9] 




phot: [10] 




phot: \20\ 


5533 


rot(Ha+[NII]): fll] 


7331 


rot(Ha+[NII]): [21] 




rot (HI): [9] 




rot(HI): [11, i9] 




density (HI): [22] 




density (HI): [8] 




phot: [11, 20] 




phot: [10] 




-25 25 

radial distance5,arcsec 




-25 25 

radfal di5tances,arcsec 



FIG. 3: Radial dependencies of the mean age and mean metallicity of the stellar populations in NGC 5533, 
which we have obtained with the long-slit spectral data from the reducer SCORPIO of the Russian 6m 
telescope. 



radial distributions of the mean age and metallicity of the stellar populations in their disks by 
analysing long-slit spectral data from the reducer SCORPIO of the Russian 6-meter telescope [l^. 
So for these galaxies we are able to use the appropriate mass-to-light ratio values obtained in the 

evolutionary synthesis models by [l^ (and re-calculated by us for the reduced Salpeter IMF so 

f j 

slightly increased with respect to the models by [14:] which used the Kroupa IMF). For the other 
two galaxies, NGC 2841 and NGC 7331, the disk stellar population parameters are unknown; for 
them we have used the empirical calibrations of the mass-to-light ratio values versus broad-band 
colours from 15| - one versus q — i in NGC 2841 and the other versus B — R in NGC 7331 (the 
broad-band colour profiles for them are also taken from the photometric survey |1Q|). 

Earlier we have found two Jarge-scale stellar disks with different characteristics of the stellar 
populations in NGC 7217 



23 



24j : inside R w 45" a thin stellar disk with the mean stellar age of 



5 Gyr and the mean metallicity of [Z/H]= —0.2 dominates. At larger distances from the center a 
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thicker and more extended star forming stellar disk with a younger, (T) = 1 — 2 Gyr, and very 
metal-poor, [Z/H]< —0.4, stellar population is seen. The border between two disks is rather narrow 
so a jump in the mass-to- light ratio is noticeable at i? ~ 50"; the matter density profile (Fig. [2]) 
demonstrates a break at the same radius. NGC 5533 is also a multicomponent stellar system, 
however this time the limits of the subsystems do not coincide in the brightness and density profiles: 
all the subsystems are more extended as presented by the gravitating matter density distribution. 
The surface photometry allows to distinguish a bulge-dominated region within R < 20" , an inner 
stellar disk at R = 20" — 60", and an outer stellar disk of low surface brightness at R > 70" 
By applying the software ULYSS [25| to the long-slit spectral data obtained with the SCORPIO of 
the 6m telescope, we have succeeded to derive the stellar populations characteristics up to « 50", 
so only for the inner stellar disk. The results are shown in Fig. [3l One can see that at i? > 20", 
or beyond the bulge-dominated area, both the age and the metallicity demonstrate fiat radial 
dependencies. So for the inner stellar disk in NGC 5533 we can use the mean values of the stellar 
population characteristics at i? > 20", namely, (T) = 2 Gyr and [Z/H]= —0.4. We have taken 
the model mass-to-light ratio just for these stellar population parameters from jl4], and then have 
added 0.15 in the logarithm to re-calculate it with the reduced Salpeter initial mass function (IMF) 
instead of the Kroupa IMF used by [l^ (the recommendation are taken from [15] ) . 

So in Fig. S] we confront the observed surface brightness profiles with those calculated from 
the matter density profiles reconstructed from the rotation curves. If the main contributor into 
the gravitational potential of a galaxy is a stellar disk, two profiles must coincide. What have we 
discovered from the Fig. We see that the different galaxies demonstrate different contributions of 
their stellar disks into the total gravity. In NGC 2841 and NGC 7217 the inner parts of both profiles 
are in agreement; in these galaxies, indeed, within R < 9 kpc in NGC 2841 and within i? < 4 kpc 
in NGC 7217, the stellar disks are the main contributors into the total gravitational potentials. 
In the outer parts of the galaxies a 'hidden mass' reveals itself - the profiles reconstructed from 
the rotation curves go strongly above the observed surface brightness profiles. But in other two 
galaxies, NGC 5533 and NGC 7331, the observed and calculated profiles diverge everywhere in the 
galaxies. In NGC 5533 the 'hidden mass' dominates in the gravitational potential over the whole 
galaxy extension - it is just the same conclusion that we have already made in our previous work 



basing on the analysis of the rotation curve shape [ll[. NGC 7331 represents a more intriguing 
situation: the observed central brightness exceeds our expectations arisen from the low visible 
rotation velocity. We think that non-circular gas motions are significant in the center of this galaxy 
so the rotation curve derived from the single long-slit cross-section along the isophote major axis 
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FIG. 4: The comparison of the observed surface brightness profiles for our four galaxies (signs), with the 
ones based on the gravitating matter surface density profiles restored from the rotation curves for the outer 
rotation curve asymptotic taken as a constant value (lines) calculated under the assumption that all the 
gravitating matter is the stellar component of the galactic disks. 

does not reflect a circular speed which characterizes the gravitational potential. Our suspicion 
about the tri-axiality of the gravitational potential shape in NGC 7331 has also been expressed by 
us earlier |26| . 

What properties of NGC 2841 and NGC 7217 may be related to the dominance of their large- 
scale stellar disks in their total gravity? We can do a couple of suggestions. First of all, both galaxies 
are classified as isolated and are included into the list of isolated galaxies by {27]. Secondly, there is 
a noticeable difference in neutral hydrogen content between NGC 2841 and NGC 7217, on one hand, 
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100 200 300 -lOO 100 200 300 400 

R, arcsec R, arcsec 



FIG. 5: The comparison of the gravitating matter density profiles reconstructed by us from the rotation 
curves with tlic outer rotation curve asymptotic taken as a constant value with the neutral hydrogen radial 
distributions in four galaxies under consideration. 

and NGC 5533 and NGC 7331, on the other hand: the latter pair of galaxies is more gas-rich (see 
the Table H]) . In Fig. [5] we compare radial profiles of the gravitating matter density reconstructed 
by us from the rotation curves, to the neutral hydrogen surface density radial distributions (at the 
logarithmic scales) for all four galaxies. While in NGC 2841 and NGC 7217 the neutral hydrogen 
contributes less than one percent of the total surface density, in NGC 5533 and NGC 7331 its 
contribution is as high as a few percent over the whole extension of the galactic disks. Moreover, 
the total character of the gas density and matter density profiles is the same: they look like a 
quasi-exponential decreasing function which is very typical just for stellar and gaseous disks of 
galaxies. The remark that the 'hidden mass' distribution in spiral galaxies is quite similar to the 
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'visible' gas distribution has been already made earlier by e.g. The fact found by us in this 

work that the gravitating matter surface density distribution agrees better with that of the stellar 
component in the gas-poor galaxies is therefore in line with the hypothesis by [29l] that the 'hidden 
mass' in spiral galaxies is in fact very cold molecular gas which cannot be observed in accessible 
electromagnetic radiation spectral domains. 



VI. ACCOUNTING FOR THE COMPONENTS OF A SPHERICAL HALO OF THE 

GALAXY 

In all previous sections we ignored the fact that much of the matter of the galaxy, in principle, 
can be concentrated in a spherically symmetric halo, which has a mass function Mgpherei''')- Initially, 
we do not know the distribution of the spherical mass Mgph(,re{r) , this feature is a free parameter 
model. 

Due to the additivity of the potential ^ (in the Newtonian approximation), it can be divided 
into two parts (the disk and a spherical component): 

#(p, Z) = ^discip, Z) + ^sphereir) , r"^ = + Z^ ■ (23) 

Then equation ([1]) in the equatorial plane {z = 0) can rewritten in as^: 

d^dics{p,Z = ^) ^ V'^jp) - GMsphere{p)/p ^24) 

dp p 

Hence, all subsequent arguments (starting with the formula ([I])) amounts to replacing the observed 
features V{p) as: 

V{p) ^ ^V^{p)-GMsphere{p)/p (25) 

Here we are free to vary the expression GMgp/iereip)/ P in the range from zero to a value approx- 
imately equal to V'^{p) — GMdisdp)/ P, where Mdisc{p) denotes the mass which is concentrated in 
the disk to a certain radius p. You can, for example, from a certain radius (e.g. the radius of the ob- 
served stellar disk to the end of the 25-th blue isophotes), put V^lp) — GMsphereip)/ P = GM^J^^^/ p, 
where Mj^°^^ is the total mass of matter in the galactic disk. This will correspond to the termination 
of the matter density in the disk at this radius (25-th blue isophotes) . 



^ Here we took into account that for the spherical part of potential we have: ^sphere{r),r = GMsphere(r) / r'^ . 
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VII. CONCLUSIONS 

In this paper we considered a new method for reconstructing of the surface density of matter 
in a flat disk of spiral galaxies. At the beginning we suppose that all matter is concentrated in 
the disk of the galaxy and apply the method for reconstructing the surface density of matter in 4 
spiral galaxies: NGC2841, NGC7217, NGC7331 and NGC5533. We compared our results with the 
distributions of the stellar matter. This comparison allows us to obtain some conclusions about 
dark matter in these galaxies. In section IVll we discuss possible accounting for the components of 
the spherical halos. 
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Appendix A: Numerical Integration of our new method 

As described in the main text, we wish to use our new method for finding distributions of surface 
density and mass throughout a galactic disk, i.e. eqs. [T2]and[T3j Due to the observational nature of 
our input data (the velocity function), our only means for solving these equations are by numerical 
methods. This appendix describes and tests our numerical code developed for this purpose. 

The fundamental integral of our method, eq. \T2\ do not have non-integrable singularities, 
however, the integrand may be ill-defined at the integration endpoints. Because of that, we use 
the open-ended second Euler-Maclaurin integration formula as our basic integration algorithm. 
Secondly, we use increasing numerical resolutions and Romberg's method to achieve the highest 
possible accuracy (see discussion on these techniques and further references e.g. in ^^). We use 
these methods with enough iterations until the result has satisfactorily converged. 

In addition to these basic integration techniques, a few other "tricks" are needed. Our basic 
data for integration will usually be a discrete set of velocity data in the disk as a function of the 
radius, p. Obviously, eq. [I2]is a function of the velocity function, but because the velocity function 
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is expressed through the intermediate variable i?o = (see eqs. [TT]fT2]l . we will often need 

values of the velocity function not explicitely tabulated by our observational data. To solve this 
issue we need to use interpolation for velocity data in between observational data points and also 
we need to make some assumptions about the behaviour of the velocity function at points beyond 
the tabulated range of observational data. 

To obtain velocity data for points in between the observational data points, we have used simple 
linear interpolation (attempts with using higher order interpolations were made, but yielded no 
increased accuracy). Furthermore, we have made the assumptions that 1) the velocity at the 
central axis is zero and 2) that the velocity beyond the outermost observed point is a constant 
function with a value equal to that of the outermost observed point. 

Finally, we have tested our numerical implementation of the new method by doing a number of 
convergence tests. We define an average error measure as : 

-.^17 ?■ I 

^^(-) = -E^^^4^- (Al) 

where x denotes the function we are investigating, n denotes the number of data points (in p) and 
N denotes the number of iterations we use in Romberg's method'*. The data in figure [6] uses this 
error measure for the two basic basic functions U{p) and (t{p) (c.f. eq. [T2]l for varying number of 
iterations. These two functions are the most basic functions of our new method and if the method 
is correctly numerically implemented, the average error must converge to zero as we increase the 
numerical resolution (number of Romberg iterations). 

As the input data for this figure we have used the test velocity function eq. [171 but it is noted 
that similar results are obtained if actual observational velocity data is used. From the figure 
it is clearly seen that both U{p) and (y{p) are converging towards zero with incresing number of 
iterations. This is exactly what we expect and we thus conclude that the numerical implementation 
of our new method is correct. 

Finally it is noted that for all numerical results presented in this paper, convergence tests have 
been made to confirm that the presented data has converged to its limiting value. 
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